Novel indicator geostatistics for water table mapping that incorporate elevation, land use, stream network and physical constraints to provide probabilistic estimation of heads and fluxes

Mapping of groundwater level observations often makes very little use of auxiliary data and is often undertaken simply by manual interpolation or ordinary kriging of the heads. Recently, a number of geostatistical methods have emerged that significantly improve estimates by incorporating the land surface elevation and groundwater flow or drawdown equations. However, at the regional scale heads are influenced by numerous other factors that cannot be considered by these methods. Such factors include the land cover type, aquifer basement elevation and upper limits to the heads (such as the land surface). Furthermore, all existing methods fail to include observation uncertainty; produce poor measures of prediction uncertainty; and assume the random field to be multi-Gaussian; that is, the spatial correlation in heads are independent of the head magnitude. To overcome these limitations and to make better use of the observation data, this paper presents a novel indicator geostatistical simulation method for mapping unconfined heads. The simulation method produces many equally probable maps and by post-processing produces quantitative uncertainty maps. Other post-processing could produce new products such as the probability of a stream having a gaining or losing hydraulic gradient and, if multiple time points are mapped, probabilistic changes in storage. To demonstrate the methodology, this paper presents an application for the Broken catchment, Victoria. The method combines a multi-variate version of kriging with external drift (KED) and a modified MarkovBayes indicator simulation algorithm to facilitate inclusion of physical constrains to groundwater head and soft data such as landuse. The KED facilitates inclusion of continuous variables that are linearly correlated with head and is used to produce a surface that results from these variables alone. As the difference between this surface and the observations were found to be spatially correlated, and approximately first and second order stationary, the head estimate was able to be refined using indicator kriging (IK) simulations. While IK was essential for inclusion of the land class data and the groundwater head constraints, it also allowed the spatial correlation to vary with the magnitude of the heads. In effect this means it relaxes an assumption required for multi-Gaussian methods such as sequential Gaussian simulations. The entire methodology was implemented within the R statistics package using the Gstat library and modified GSLib algorithms. The source code will be made publicly available with a forthcoming journal paper. The study area comprised of both the Broken River and Broken Creek catchments with a 20 kilometre buffer to minimise boundary artefacts. Groundwater observations comprised of data from The Department of Primary Industries, Victoria; The Department of Sustainability and Environment, Victoria; and Department of Water and Energy, NSW. Importantly, the land class was found to be statistically important in estimating heads and the spatial correlation was not multi-Gaussian, thus challenging the use of standard kriging methods. All simulation head maps honoured all of the constraints. However, because of large nugget values and observations not located at the centre of grid cells, grid cells estimates can differ from observed water levels. Despite this implementation issue, overall the methodology produced highly plausible water table maps that increase the information extracted from observation bores by use of a very wide range of quantitative and qualitative spatial data.